package org.darkbits.cyklon.util;

public class Transform {
	public static Complex[] fft(Complex in[]) {
		if (!isPowerOfTwo(in.length)) {
			throw new IllegalArgumentException("Length of argument for FFT is not a power of two");
		}
		
		return recursiveFft(in, 1.0f);
	}
	
	public static Complex[] fft(float in[]) {
		if (!isPowerOfTwo(in.length)) {
			throw new IllegalArgumentException("Length of argument for FFT is not a power of two");
		}
		
		Complex cin[] = new Complex[in.length];
		
		for (int i = 0; i < cin.length; i++) {
			cin[i] = new Complex(in[i], 0.0f);
		}		
		
		return recursiveFft(cin, 1.0f);
	}
	
	public static Complex[] ifft(Complex in[]) {
		if (!isPowerOfTwo(in.length)) {
			throw new IllegalArgumentException("Length of argument for IFFT is not a power of two");
		}
		
		Complex out[] = recursiveFft(in, -1.0f);
		
		float scaling = 1.0f / in.length;
		for (int i = 0; i < out.length; i++) {
			out[i] = out[i].mult(scaling);
		}
		
		return out;
	}
		
	private static Complex[] recursiveFft(Complex a[], float sign) {
		int n = a.length;
		int n2 = n / 2;
		if (n == 1) {
			return a;
		}
		Complex on = new Complex((float)(sign * 2.0 * Math.PI / n));
		Complex o = new Complex(1.0f, 0.0f);
		Complex a0[] = new Complex[n2];
		Complex a1[] = new Complex[n2];
		for (int i = 0; i < n2; i++) {
			a0[i] = a[i * 2];
			a1[i] = a[i * 2 + 1];
		}
		Complex y0[] = recursiveFft(a0, sign);
		Complex y1[] = recursiveFft(a1, sign);
		Complex y[] = new Complex[n];
		for (int k = 0; k < n2; k++) {
			y[k] = y0[k].add(o.mult(y1[k]));
			y[k + n2] = y0[k].sub(o.mult(y1[k]));
			o = o.mult(on);
		}
		return y;
	}
	
	private static boolean isPowerOfTwo(int i) {
		return i > 0 && i == (i & -i);
	}
	
	public static float[] hanning(int size) {
		float r[] = new float[size];
		
		if (size == 1) {
			r[0] = 1;
		} else {		
			for (int i = 0; i < size; i++) {
				double n = i / (size - 1.0);
				r[i] = (float)(0.5 + 0.5 * Math.sin(n * Math.PI));
			}
		}
		
		return r;
	}
	
	public static float[] powerSpectrum(float in[]) {
		Complex cout[] = fft(in);
		
		float result[] = new float[in.length / 2 + 1];
		
		float scaling = 2.0f / (in.length * in.length); 
		for (int i = 0; i < result.length; i++) {
			result[i] = cout[i].abs2() * scaling;
		}
		result[0] /= 2.0f;
		result[result.length - 1] /= 2.0f;
		
		return result;
	}
}
